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A PROGRAMMED MATHEMATICAL MODEL 
TO SIMULATE THE BENDING OF RADIO 
WAVES IN ATMOSPHERIC PROPAGATION 


B. Rosenbaum 
N. Snow 


ABSTRACT 


The ray path of an RF signal undergoes bending in the interval of atmos- 
pheric propagation between a satellite and an earth-based tracking station. 
Resultant refraction is a perturbation on observed elevation angles and Doppler 
shifted signals. As a practical tool for the computation of corrections a mathe- 
matical model simulating the tracing of the refracted ray path has been pro- 
grammed for use by the Goddard Space Flight Center computer facilities. 
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A PROGRAMMED MATHEMATICAL MODEL 


TO SIMULATE THE BENDING OF RADIO 
WAVES IN ATMOSPHERIC PROPAGATION 


I. INTRODUCTION 

The ray path of an RF signal undergoes bending in the interval of atmos- 
pheric propagation from a satellite to an earth -based tracking station. The 
refraction introduces a perturbation on tracking measurement of elevation angle 
and measurement of Doppler shift and its conversion to range rate. Correction 
of the effect can be performed through a calculation of the finite refraction 
angles between the line-of -sight and the ray path at both the emitter and receiver. 

This document presents the mathematical model and its associated computer 
program to simulate the ray tracing of a signal path in a spherically stratified 
atmosphere. While many such routines performing similar traces exist (Refer- 
ence 1), they are not adapted to the newly installed IBM 360 computers at Goddard 
Space Flight Center (GSFC). This program has been checked out and run at GSFC 
on computer models IBM 360/91 and IBM 360/95. An exponential model is used 
for the tropospheric refractivity profile and the Chapman function is used for the 
ionospheric profile. The refraction effect on elevation angle attributable to the 
troposphere and ionosphere can thus be calculated from a single routine. The 
computation procedure is essentially that of the National Bureau of Standards 
tropospheric ray tracing program (Reference 2). A sample calculation is in- 
cluded in Chapter IV, Users Guide. 


n. REFRACTIVITY MODEL 

The following equations were used in the construction of the refractivity 
model. 

1. Relation of Refractivity to the Index of Refraction 

The relationship of refractivity, N, to the index of refraction, n, is expressed 
by the equation: 


N = (n - 1) x 10 6 
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2. Tropospheric Refractivity Profile 


The exponential tropospheric refractivity profile is expressed as: 
N(r) = N(r 0 + h) = N 0 e“ ch 

where 

0 < h < SO km 


3. Ionospheric Refractivity Profile — Chapman Model 

The refractivity in the ionosphere is a function of free electron density, N e 
and signal frequency, f , as expressed in the following equation: 

-40.38 X lo 6 N e (r 0 + h) 


where 

h > 50 km 


The altitude profile for electron density is given by the Chapman function: 
N e = N m exp 1/2(1 - z - e“ z ) 

where 

■ ( h - h J 


m. METHOD OF COMPUTATION 

The ray tracing procedure divides the atmosphere into a series of thin 
laminations (Figure 1). The refraction between consecutive layers is then fol- 
lowed. Snell's law, which for a spherically stratified medium is given by: 



RAY PATH 



Figure 1. The Geometry of-the Ray-Tracing 
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is the basis for the calculation. This is rewritten for computation purpose 
as: 


sin 2 ( 8.J1\ ) 



2 sin 2 


(V 2 ) 



(N 0 - Nj) 


n. 


x 10' 


r6 



The incremental bending angle, At, is determined by the general formula: 

An 

At — cos 8 

n 

Between the i th and (i + l) th lamination the bending, At. , is calculated by 
the expression: » 


Ar i " 


(Ni - N i+1 ) , \ 


(n 


The total bending is given by: 


-v = 

i = l 

The error in the computation of r L is made small by making the laminations 
sufficiently thin. 

The refraction angle is calculated from the formula: 

cos t. - sin t. % tan 6 . ~ n./n 0 
tan e A - ( n -/ no ) tan 0 O - sin r L - cos t. tan 0 i 

and the refraction angle (h can be calculated by the formula: 
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tan S- 
1 


n 0 /n i - cos r 1 - sin t 1 tan 8 0 
sin r. - cos r x tan $ 0 + (n 0 /n A ) tan 9 1 


The relation, - e ± + S l9 can be used as a check on the computation. 

IV. USERS GUIDE 

1. Introduction 

This program computes the bending angle of radio waves in atmospheric 
propagation. 

2. Programming System 

The program is written in double precision, FORTRAN IV compiler language. 
It is operational on the IBM 360/Model 91 and Model 95 computers. 

3. Required Input/Output Units 

No special equipment configurations are required. Input is a source card 
deck and data cards. Output is a standard printout. 

4. Restrictions 

There are no restrictions on this program. 

5. Input 

The following values are read into the computer in a DIO. 2 format. 

(1) The reciprocal scale height of the troposphere (c), 

(2) The maximum electron density (N m ), 

(3) The initial value of the refractivity (N 0 ), 

(4) The signal frequency in megahertz (f), 

(5) The scale height of the ionosphere in kilometers (2H), 

(6) The height of N m , (hj, and 
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(7) The initial value of the angle between the ray path and the* radius vector, 
in radians (& 0 ). 

A sample input data card is shown below. 


15.6SB-G2 


(0 


COCO 


CO 


34.45B+G1 


(3) 


1 . 3<ED+t?Sj 


(4) 


(.5) 


Ua) 


OOf .401+0^ 


(7) 


000000. 0 0 op . 0 0 0,0 0 0 0 0 p 0 0 0 0 0 

1 2 3 4 5 6 7 8 3 10 pi 12 13 14 15 IE 17 18 19 20 21 22 23 24 25 26 27 23 29 3(7 

i i 1 1 1 1 1 in 1 1 1 n ii 1 1 n 1 1 1 1 1 1 i 


0 0 0 0 0 0 0 0 010 0 0 0 0 0 0 0 I jo 1 I 0 0 0 0 0 

31 32 33 34 3 5 36 37 38 39 40 41 42 43 44 45 4 6 47 48 43 50 31 52 S3 54 5 5 56 57 58 59 60 

11 1 1 1 1 11 1 1 11 1 11111111111111 


00000 00 


oooooooooo 


SI £2 63 64 65 6 S 67 58 69 70f71 72 73 74 73 76 77 7 6 79 M 

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 11 1.1 1 


The following Job Control Language cards are required to successfully 
execute the program. 

C.C.* 

1 

//DSNAS033bJOBb(G7001H1310,T,DA0023, 001005), 608, MSGLEVELKL 

//bEXECbFORTRAN 

//SOURCE. SYSINbDDb* 

SOURCE LANGUAGE PROGRAM 


C.C.* 

1 

/* 

/ /bEXECbLINKGO 
//GO.DATA5bDDb* 


DATA CARDS 


C.C * 

1 

/* 

* C.C. = Card Column 
b = blank column 
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6. Output 


The values of the following parameters are printed as the input value of H 
is incremented by the program to its predetermined maximum (3000 km). For- 
mat is indicated following the parameter. 

H (altitude in kilometers) F 7.1 

N (r ef ra ctivity ) E 1 5 . 5 

N. “ N i + i (difference between successive values of N) E15.5 

8 (angle between ray path and radius vector, in 

milliradians ) E 1 5 . 5 

Ar (discrete portion of the bending angle, in milliradians) E15.5 

r (bending angle, in milliradians) E15.5 

€ (refraction angle at the earth surface, in milliradians) E15.5 

S (refraction angle at the satellite, in milliradians) E15.5 

e + S (the sum of epsilon and delta, in milliradians) E15.5 

In addition, initial values and constants used in the equations are printed. A 
sample output follows. 
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RENDING DF RADUl KAVFS IN ATRnSPHF» I C RRfiRAGAT TON 
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7. Program Listing 

A complete listing of the program follows: 

G LEVEL 1, MOD l MAIN DATE = 6812* 19/59/1? 

DOUBLE PRECISION THETA! 125 ) ,DELTAU< 1 25) .EPS! I25).SN!125).AN(1?S), 
1HAUXC6) ,TAU( 125) , GAMMA! 1 25 ), H ( 1 25) 

DOUBLg PRECISION OARS IN . DEXP , DSQRT, D ATAN . DTAN , DCOS . THET AO 
OOUBLE PRECISION COTHE.TATHE 

DOUBLE PRECISION D , E , SUM , C , G AMI ,E PA , A , 8 , DE N l , OEN2 , DEN 3 , ANUM 1 , 

1 ANUM2. ANUM3 , SNO 

DOUBCE PRECISION ANM , HM , HH , ANO . F , RO . ANE » Z . EX 1 . EX 2 • APG » D I ff 
EQUIVALENCE ( HAUX ( 1) . H ( 1 ) ) 

DATA HAUX/O.DO, . 1 00 , . 300 , . 5DC , . 700 , I .000/ 

100 FORMAT (8010.2) 

DO 40 1=7,21 
40 H( I > = H( 1-1 >+.2D0 
DO 42 1=22,33 
42 H ( I )=H( 1-1 J+.5D0 
DO 43 1=34,38 

43 H{ I)=H( I— 1 ) + 1.0 
H( 39 )=20 .ODO 

DO 44 1=40,108 

44 H( I )=H( 1-1 ) +10. ODO 
DO 45 1=109,122 

45 HCI)=H( 1-1 1+20.0DO 
H ( 1 23 ) = 1 000 .ODO 

H ( 1 241=2000.000 
H ( 1251=3000. ODO 

C PROGRAM STOPS WHEN THFRE ARE NO MORE DATA CARDS TO BE READ IN 
200 RE AD ( 5, 1 00 , END = 1 22) C , ANM , ANO , F , HH , HM , THET AO 
ANM=4. 248550 + 1 1 
TATHE=DTAN( THET AO) 

CQTH£=DCOS( THET AO ) 

T 0=THET A 0*1 .0+3 

R0=6373.01 5D0 

SN0=1 .D0+AN051 .D— 6 

EX1=2.D0*(DSINCTHETA0/2.D0 ) ) **2 

DO 62 1=1,125 

AH=H ( I ) 

IF ( AH.GT .50 .0 ) GO TO 53 
ANC I )=ANO*DEXP(-C*H{ I ) ) 

SN (I ) = 1 . D0 + AN (I)*l.D— 6 
GO TO 60 

58 Z=(Ht I)-HM)/HH 

ANE=ANM*DEXP( ( 1 . DO- Z-DEXP ( -Z ) 1/2.D0 ) 

AN (I ) = <-40.38D0*ANE*l .D+6)/F**2 
SN ( I ) = 1 .DO + AN!I)*l.D-6 
60 CONTINUE 

EX2=( ANO-ANI T ) ) /SN(t) *1.D-6*C0THE 

ARG=DSQRT(R0/(2. + (R0+H( I ) ) ) * ( E X 1 +H(I)/RO -EX2)) 

THETA { I)=2.D0*DARSIN(ARG) 

62 CONTINUE 

WRITE <6 .105) ANM , HM »HH,AN0»T0»C»F,R9 
WR I TE ( 6 , 108) 

WRITEfB, 110) 

I LCNT= 1 4 
THETA! 1 )=THETAO 
EPS< 1 )=0.0D0 
GAMMA! 1)=0. ODO 
TAU! 1 )=0 . ODO 
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